
function [ys,rw,res]=marray_lowess(x,y,f,nsteps,delta)
%[ys,rw,res]=lowess(x,y,f,nsteps,delta)
% x: must be ordered from smallest to largest. 
% y: cordinates of the points on the scatterplot.
% f: is the fraction of points used to compute each fitted value
% as f increases the smoothed values become smoother, choosing f from 0.2 to
% 0.5 should be solve most of the problems
% nsteps: the number of iterations in the robust fit. nsteps=0 nonrobust fit is 
% returned. nsteps=2 should server most purpose.
% delta: nonnegative parameter which may be used to save computations; if n is less than
% 100, set delta to 0, if n greater than 100 you should try by yourself, for
% example, delta=length(x)/k, k=50.
%
n=length(y);
if n<2
   ys(1)=y(1);
end
ns=max(min(fix(f*n),n),2);
for iter=1:nsteps+1
   nleft=1;
   nright=ns;
   last=0;
   i=1;
  while last<n
   while nright <n
      d1=x(i)-x(nleft);
      d2=x(nright+1)-x(i);
      if d1>d2
        nleft=nleft+1;
        nright=nright+1; 
      else
        break;
      end  
   end
   
   
   if iter>1
    [ys(i),res,ok]=marray_lowest(x,y,n,x(i),nleft,nright,iter>1,rw);
   else
     [ys(i),res,ok]=marray_lowest(x,y,n,x(i),nleft,nright,iter>1);
   end
  if ~ok
      ys(i)=y(i);
   end
   if last<i-1
      denom=x(i)-x(last);
      for j=last+1:i-1
         alpha=(x(j)-x(last))/denom;
         ys(j)=alpha*ys(i)+(1-alpha)*ys(last);
      end
   end
   last=i;
   cut=x(last)+delta;
   for i=last+1:n
      if x(i)>cut 
         break;
      end
      if x(i)==x(last)
         ys(i)=ys(last);
         last=i;
      end
   end
   i=max(last+1,i-1);
end

for i=1:n
   res(i)=y(i)-ys(i);
end

if iter>nsteps
   break;
end

for i=1:n
   rw(i)=abs(res(i));
end
rw=sort(rw); %Check

m1=n/2+1;
m2=n-m1+1;
m1=ceil(m1);
m2=ceil(m2);
cmad=3*(rw(m1)+rw(m2));
c9=0.999*cmad;
c1=0.001*cmad;
for i=1:n
   r=abs(res(i));
   if r<=c1
      rw(i)=1;
   elseif r>c9
      rw(i)=0;
   else
      rw(i)=(1-(r/cmad)^2)^2;
   end
 end
end

      

   

   
         
   


